Stress Relaxation in Entangled Polymer Melts 
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We present an extensive set of simulation results for the stress relaxation in equilibrium and step- 
strained bead-spring polymer melts. The data allow us to explore the chain dynamics and the shear 
relaxation modulus, G(t), into the plateau regime for chains with Z = 40 entanglements and into 
the terminal relaxation regime for Z = 10. Using the known (Rouse) mobility of unentangled chains 
and the melt entanglement length determined via the primitive path analysis of the microscopic 
topological state of our systems, we have performed parameter-free tests of several different tube 
models. We find excellent agreement for the Likhtman-McLeish theory using the double reptation 
approximation for constraint release, if we remove the contribution of high-frequency modes to 
contour length fluctuations of the primitive chain. 

PACS numbers: 83.80.Sg (Polymer melts), 83.10.Rs (MD simulation rhcology), 61. 25. he (liquid structure 
polymer melt) 



High molecular weight polymeric liquids display re- 
markable viscoelastic properties [HQ ■ Contrary to glassy 
systems, their macroscopic relaxation times are not due 
to slow dynamics on the monomer scale, but arise from 
the chain connectivity and the restriction that the back- 
bones of polymer chains cannot cross while undergoing 
Brownian motion. Modern theories of polymer dynam- 
ics 0, 0| describe the universal aspects of the viscoelastic 
behavior based on the idea that molecular entanglements 
confine individual polymers to tube-like regions in space 

i| . Forty years of research have led to a complex relax- 
ation scenario based on a combination of local Rouse dy- 
namics, reptation, contour length fluctuations, and con- 
straint release Q. The development and validation of 
a quantitative, microscopic theory crucially depends on 
the availability of experimental and simulation data for 
model systems. 

Entangledpolymers are studied experimentally using 
rheology [3, 0, , dielectric spectroscopy [§] , small-angle 
neutron scattering 0, [l(|, and nuclear mag netic reso- 
nance [III [lH . Computer simulations [l3l4lrj | offer some 
advantages in the preparation of well-defined model sys- 
tems and the simultaneous access to macroscopic behav- 
ior and microscopic structure and dynamics. In par- 
ticular, the recently developed primitive path analysis 
(ppa) [n- l24j reveals the experimentally inaccessible 
mesoscopic structures and relaxation processes described 
by the tube model and allows parameter-free comparisons 
between theoretical predictions and data. However, the 
long relaxation times pose a particular challenge to com- 
putational techniques. Here we present simulation re- 
sults for model polymer melts in equilibrium and after a 
rapid, volume-conserving uni-axial elongation, where we 
have been able to follow the full relaxation dynamics deep 
into the entangled regime. The data allow us to perform 



the first parameter-free test of the predictions of tube 
models for dynamical properties, to pinpoint a problem 
in the current theoretical description, and to validate a 
suitable modification. 

Our numerical results are based on extensive molecu- 
lar dynamics (MD) simulations of bead-spring polymer 



melts |13(. Each chain is represented as a sequence of 



beads connected by finite-extensible, non-linear (FENE) 
springs and interacting via the repulsive part of the 
Lennard- Jones 12-6 potential (LJ). The energy scale is 
set by the strength of the LJ interaction, e, while the dis- 
tance scale is set by the monomer size, a. The basic unit 
of time is r = o^m/e) 1 / 2 , where m is the mass of each 
monomer. The equations of motion are integrated using 
the LAMMPS MD simulation package with a veloc- 
ity Verlet algorithm and a time step St = 0.012r. The 
temperature, T = e/fcs, was kept constant by weakly 
coupling the motion of each bead to a heat bath with a 
local friction T = 0.5t _1 . 

We have studied seven entangled polymer melts of M 
chains of N beads with M X N = 5000 x 50, 2500 x 100, 
400 x 175, 200 x 350, 200 x 700, 400 x 1000, and 320 x 3500 
each at a monomer density p = 0.85er -3 . Using the most 
refined PPA estimate of the rheological entanglement 
length for this model of N e = 85 ± 7 [36j , the investigated 
systems span the range from unentangled (Z = N/N e w 
0.6) to highly entangled (Z = N/N e w 41). The Rouse 
time was previously determined as tr = 1.5r./V 2 [l3j|, en- 
tangemcnts effects become relevant around r e = Tu(N e ), 
and the maximal relaxation times of entangled systems 
are expected to be on the order of r° = r e Z 3 = trZ . 

The melts were generated and equilibrated following 
the procedure outlined in Auhl et al. [l6j . Techni- 
cally, the largest challenge is the reliable extraction of the 
macroscopic, viscoelastic behavior [26j, |27[. Data were 
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FIG. 1: (color online) (a) Normal tensions, a(X,t) for step- 
strained melts N = 50 (black), 100 (blue), 175 (green), 350 
(red), 700 (cyan), 1000 (violet), and 3500 (orange) and elon- 
gation A = 2 (□), 3 (o), and 4 (A), (b-d) Green- Kubo shear 
relaxation moduli, G(t), (this work: large solid O with solid 
line; Ref. [2(|: no symbol, solid line) compared to extrapola- 
tions, G(\,t) = a(\,t)/h(X,t), from the non- linear response 
using the (b) classical, (c) Doi-Edwards, (d) slip-tube damp- 
ing function (same symbols as in (a)). Colored ticks indicate 
the Rouse time of the corresponding systems. 

recorded in equilibrium as well as out of equilibrium after 
a step-strain. Strained melts were prepared by subject- 
ing equilibrated initial conformations to rapid (T^e/ G 
[120t, 36000r]), volume-conserving (X x X y X z = 1), elon- 
gational deformations with (X x = A, X y =X Z =l/\/A) for 
A ranging from 1.5 to 4.0 well outside the linear elastic 
regime. Deformations in this range are typical for many 
applications of polymeric systems and large enough to 
generate a measurable elastic response for the present 
system sizes. In the ideal case, strain should be intro- 
duced instantaneously. To check the dependance of our 
results on T^e/i we have varied the deformation time in 
one case ( N = 700, A = 3.0) by a factor of 200. In the 
following, for deformed systems, t = is fixed to the mid- 
dle of the deformation period, i.e. data are recorded for 
t > Td e f/2. To reduce finite- Td e f artifacts, we typically 
discard data from the initial 3Td e f- 

The longest simulations were run up to 2 x 10 9 time 
steps and sufficient to reach the plateau regime for 
our longest chains and to completely relax the others 
(2 x 10 9 x 0.012t ps t r {N = 4000) ps r d {Z = 12.5)). The 
total numerical effort corresponds to about 5 million sin- 
gle core CPU hours. We recorded block-averages of the 
microscopic stress tensor u a p{t) = Fij,a.fij,p) /V 

at intervals of 1.2r. The latter sum is over all pairs i,j 
of interacting beads, a, f3 are Cartesian indices, and F, r 
and V denotes force, separation and volume, respectively. 
Furthermore, we stored melt conformations at intervals 
of 120r for further analysis of the chain conformations. 



Results for the relaxation of the normal tension a(t) = 
&xx — \ { G yy + a zz) are presented in Fig. QJi. We observe 
a clear non-Newtonian behavior with a stress relaxation 
extending over many orders of magnitude in time after 
the end of the deformation period of the sample. As 
expected, there is a strong increase of the terminal relax- 
ation time with chain length and the gradual formation 
of an intermediate plateau in the stress relaxation for the 
longest chains studied. The maximal relaxation time is 
independent of the total deformation. 

Figure HJb-d) show comparisons between 
a(X,t)/h(X,t) for the step-strained melts to the 
linear shear relaxation moduli G(t) obtained by Likht- 
man et al. and ourselves via the Green-Kubo 

relation G(t) = V {a a p{t)a a p(Q)) / k B T with a ^ /? 
the from stress fluctuations in unstrained, equilibrated 
melts. Available expressions for damping functions 
h(X) are refinements of the stress-strain relation 
h(X) — A 2 — A" 1 predicted by classical rubber elas- 
ticity theory (Fig. [Tb). The Doi and Edwards [HI 
damping function (Fig. [T]i) includes the dynamics of 
a uniform chain retraction inside the stretched tube. 
The Rubinstein-Panyukov slip-tube damping function 
0, h{X) = (A 2 - A -1 ) / (0.74A + 0.61A -1 / 2 — 0.35), 
accounts for non-affine tube deformations and the 
asymptotic chain length redistribution inside the tube 
(Fig. [Hi). For times t < tr(N), the presence of 
additional relaxation processes prevents a systematic 
extraction of G(t) from normal tensions measured in the 
non-linear regime. Empirically, the time-independent 
slip-tube expression works surprisingly well. For times 
t > Tfl(iV), the differences between the slip-tube and the 
Doi-Edwards damping functions are small. Both result 
in a satisfactory data collapse and good agreement with 
the Green-Kubo results, suggesting that they capture 
the non-linear effects at large strains with reasonable 
accuracy. The step-strain data included in Figs. [2] and [3] 
correspond to times t > 0.25tr(N). 

In Figure [5Jt we show a comparison of the simulation 
results for G(t) to the Rouse model predictions 0, 
for unentangled systems. Our results confirm the ex- 
pectation that the Rouse model quantitatively describes 
the chain length independent early-time stress relaxation 
with G(t) ex t~ x / 2 as well as terminal stress relaxation 
in systems where the chains are too short to be entan- 
gled. For longer chains, entanglements start to affect the 
behavior beyond a material-specific, characteristic time 
r e ps 10 4 t with a gradual formation a plateau in the 
stress relaxation reached by our longest chain systems 
with Z = 41. For the terminal stress relaxation of sys- 
tems with Z = 0(10) we have reliable data extending 
about one order of magnitude below the plateau level. 
This is sufficient to allow for a meaningful comparison to 
current theories. In particular, we are not restricted to 
comparing the ability of different theories to fit the data. 
Rather, we can carry out absolute, parameter-free com- 
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FIG. 2: (color online) Comparison of the measured relaxation 
moduli to the predictions of various theories. Symbols and 
colors for the simulation data as in Fig. [JJ with Green-Kubo 
data shown as large solid O an d d. Theoretical predictions 
are shown as thick lines using the same color code. Thin lines 
indicate the uncertainty in the theoretical predictions due to 
the uncertainty of the PPA entanglement length N e — 85 ± 7. 



parisons using the result N e = 85 ± 7 [24| of the primitive 
path analysis and the known Rouse friction of the model. 

Likhtman and McLeish (LM) [3l| assembled the effects 
of (i) early-time Rouse relaxation, (ii) tension equilibra- 
tion along the contour of the primitive chains, (iii) repta- 
tion, (iv) contour length fluctuations, and (v) constraint 
release into a closed functional form, 
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where fi(t) and R(t) account for single- and multi-chain 
relaxation processes of the tube model. In their absence 
(n(t) = R(t) = 1), the formula describes a crossover from 
the early time Rouse relaxation G(t) oc t -1 / 2 to a plateau 
G% = | £ 7^. The key quantity of the tube model is the 
single-chain memory function, [i(t), for the fraction of the 
primitive chain which has not escaped from its original 
tube after a time t. Comparisons to the data neglect- 
ing constraint release (R(t) = 1) are shown for the orig- 
inal Doi-Edwards model accounting only for reptation 
(Fig. [2b), Doi's approximate inclusion of the effect of 
contour length fluctuations combining reptation dynam- 
ics with the maximal relaxation time from [3l"j (Fig. f2t) , 
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FIG. 3: (color online) Comparison of the measured relaxation 
moduli to Eq. (jTJ) using (a,b) the independently measured au- 
tocorrelation function of the primitive chain end-to-end vec- 
tors to estimate of the tube memory function, fj,(t); (c,d) 
our proposition Eqs. 121 II) for removing high-frequency modes 
from the Likhtman-McLeish theory [3l[ of contour length fluc- 
tuations with a = 1.7. Symbols and colors as in Fig. [2] 



and the full LM theory [31( of contour length fluctuations 
and reptation (Fig. [2^). For the comparisons in Figs. [2ji 
and f we have included the effect of constraint release in 
the double reptation [32| approximation R(t) = yu(i). 

The overall agreement between our data and the more 
advanced versions of the tube model is fairly good. In- 
terestingly, the more sophisticated LM theory seems to 
work less well than Doi's approximation when combined 
with the LM estimate of the maximal relaxation time. 
Yet, from our rhcological data alone, it is hard to clearly 
identify the relevance and the quality of the theoreti- 
cal description of the various relaxation processes. For 
example, one might (as we believe, erroneously; sec be- 
low) conclude, that the double reptation approximation 
strongly overestimates the contribution of constraint re- 
lease to the stress relaxation (Figs. [2fe and f) or that 
constraint release is inefficient for Z < 4 (Figs. [2]: and 
d). Obviously, fitting the various theories to the data 
would only obscure their shortcomings. 

To draw definite conclusion on how to improve the the- 
ories, we discriminate between three possible sources of 
error: (i) the functional form of Eq. (|T|), (ii) the treatment 
of reptation and contour length fluctuations underlying 
the single-chain memory function and (iii) the treat- 
ment of the multi-chain effect of constraint release via the 
double reptation approximation, R(t) = fi(t). For long 
chains and under the assumption that the escaped chain 
sections equilibrate completely, fi(t) equals the autocor- 
relation function of the chain end-to-end vectors Q , For 
shorter chains, it is more suitable to consider the end- 
point motion of the primitive chains, defined as the aver- 
age of the chain conformation over a period of r e [33| . 
This correlation function is easily accessible from our 
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equilibrium simulations and is not affected by constraint 
release Q • The comparison between the measured relax- 
ation moduli and those predicted from Eq. ([IJ using the 
measured fi(t) together with R(t) = 1 and R(t) = fi(t) is 
shown in Figs. [3] a and b respectively. For the full the- 
ory the agreement is excellent, supporting the utility of 
both the Likhtman-McLeish functional form of the shear 
relaxation modulus and of the double reptation approxi- 
mation for constraint release. 

The shortcomings of the LM description apparent in 
Fig. |2J and in the rheological study by Liu et al. Q 
must thus be related to the central part of their theory, 
the estimation of the time dependence of fi(t) under the 
combined influence of reptation and contour length fluc- 
tuations. A possible explanation is a double-counting 
of the effect of short-wavelength (p > Z) modes in the 
Rouse relaxation part of Eq. ([1]) and in fi(t). LM extrap- 
olated /i(t) to the continuum limit, resulting in a decay 
on time scales t < r e , where the motion of the primitive 
chain should be negligible. To correct for this, we have 
removed from the CLF part of fi(t) the contribution of 
modes with a relaxation time shorter than a 4 r e : 



ment length was determined independently via a topo- 
logical analysis [I?], 24 1- We find excellent agreement 
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where p* = \J Z/10 and e* are defined as in the orig- 
inal LM theory, which is recovered in the a = limit 
of the above expressions. Fig. [3ji shows that we ob- 
tain significantly improved agreement between theory 
and our data for values a = 0(1). Interestingly, 
this corresponds to a constant offset of Z = Zlm — 
■A7(«*t.) Zr °> 3 ° £ 6 5/ 4 exp(-et)de = Z LM - 1.22a. This view 
is in qualitative agreement with arguments put forward 
by van Ruymbckc et al. [34J to consider, within the orig- 



inal LM theory, chains with virtual extensions of length 
N e resulting in an increase the relaxation time of the 
outermost "real" chain segment to r e . 

To summarize, we have presented an extensive set of 
simulation results for the equilibrium and relaxation dy- 
namics of entangled model polymer melts. In particu- 
lar, we explored G{t) into the plateau regime for chains 
with Z = 41 and into the terminal relaxation regime for 
Z < 10 and compared our data to predictions of dif- 
ferent versions of the tube model. These comparisons 
did not involve any free parameters, since the cntangle- 



for the Liktman-McLcish theory using a corrected tube 
memory function and the double reptation approxima- 
tion for constraint release, demonstrating that the prim- 
itive path analysis of the microscopic structure endows 
the tube model with predictive power for d yna mical pro- 
cesses. The use of more elaborate schemes [351 ] for treat- 
ing constraint release and predicting the function R [jti(i)] 
should lead to even better agreement. 
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